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Abstract 

A technique for computing an ILU preconditioner based on the FAPINV algorithm is 
presented. We show that this algorithm is well-defined for H-matrices. Moreover, when used 
in conjunction with Krylov-subspace-based iterative solvers such as the GMRES algorithm, 
results in reliable solvers. Numerical experiments on some test matrices are given to show 
the efficiency of the new ILU preconditioner. 
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1. Introduction 



Consider the linear system of equations 

Ax = b, (1) 

where the coefficient matrix A G R nxn is nonsingular, large, sparse and x,b £ M n . Such 
linear systems are often solved by Krylov subspace methods such as the GMRES |24| and 
the BiCGSTAB [32j methods. In general, the convergence of the Krylov subspace methods 
is not guaranteed or it may be extremely slow. Hence, the original system ([1]) is transformed 
into a more tractable form. More precisely, to obtain good convergence rates, or even to 
converge, Krylov subspace methods are applied to the left preconditioned linear system 

MAx = Mb, 
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or to the right preconditioned linear system 

AMy = b, x = My, 

where the matrix M is a proper preconditioner. Two-side preconditioning is also possible 

There are two general ideas for constructing a preconditioner. The first one is to find 
a matrix M such that M = G _1 , where G approximates A in some sense. In this case, 
M should be chosen such that AM (or MA) is a good approximation of the identity 
matrix. The best-known general-purpose preconditioners in this class are those based on 
the incomplete LU (ILU) factorization of the original matrix. Let the matrix A admits 
the LDU factorization A = LDU, where L and U T are lower unitriangular matrices and 
D = diag(di, d,2, • • • , d n ) is a diagonal matrix. Here G = LDU, where L and U T are 
sparse lower unitriangular matrices which approximate L and U T , respectively, and D is a 
diagonal matrix which approximates D. The ILU preconditioners are very much effective 
in increasing the rate of convergence. Their main drawback is the possibility of breakdowns 
during incomplete factorization process, due to the occurrence of zero or small pivots (or 
appearing nonpositive pivot elements for symmetric positive definite (SPD) matrices). In 
[21] . Meijerink and van der Vorst have shown that this type of preconditioning exists for 
M-matrices. Then, Manteuffel in [22J has extended this result to the H-matrices. We recall 
that the matrix A = (ay) is an M- matrix if ay < for all a ^ j, A is nonsingular and 
A^ 1 > 0. Moreover, A is an H-matrix if its comparison matrix A = (ay) is an M-matrix 
where 

a ij | j i 7^ 3i 

O'ii | ) i — j- 

For general matrices, there are some ways to guard against appearing zero or very small 
pivots, see for example [5| 123(115]. Another drawback of the ILU preconditioners is the lack 
of inherent parallelism. Many researchers have made effort to improve the accuracy and 
the degree of parallelism of the ILU preconditioners in the past [29j [25l [26]. [27^ [28^ [3]. [4]. [20] . 

The second idea for constructing a preconditioner is to find a matrix M that directly 
approximates A^ 1 (M ~ A" 1 ). In this case, in practice we do not need to compute 
AM (or MA) explicitly, because when the Krylov subspace methods are used to solve 
a preconditioned system, only the matrix-vector product is required. One drawback of 
many sparse approximate inverse techniques is their high construction cost, unless the 
computation can be done efficiently on parallel computers. One approach in this class is 
to compute sparse approximate inverse in the factored form. Here from LDU factorization 
A = LDU, we have A" 1 = ZD^W where W = L~ l and Z = U~ l . If G = ZDW, where 
WfaW, Z » Z, and D w D _1 , in which W and Z T are sparse lower unitriangular matrices 
and D is diagonal matrix, then G may be used as a preconditioner for system ([T]) and is 
called factored sparse approximate inverse. Here W and Z are called sparse approximate 
inverse factors of A. There are several algorithms to compute a factored sparse approximate 
inverse of a matrix. Among them are the FS AI algorithm |11( [T2] , the AIB algorithm [29J , 
the AINV method [8], and the FAPINV algorithm [El Mi [13 E31 H El] . 

For SPD matrices, there exists a variant of the AINV method, denoted by SAINV (for 
Stabilized AINV), that is breakdown- free [6J. This algorithm is also presented with the 
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name of AINV-A, independently, by Kharchenko et al. in [10]. Benzi and Tuma in [5] 
have introduced an ILU factorization based on the SAINV algorithm. In the proposed 
algorithm the L factor of LDL T factorization of A can be obtained as a by-product of the 
A-orthogonalization process used in the SAINV algorithm, at no extra cost. Rezghi and 
Hosseini in [23] , have shown that a similar algorithm free from breakdown can be established 
for nonsymmetric positive definite matrices. 

The main idea of the FAPINV algorithm was introduced by Luo [15} [T6l 117] . Then the 
algorithm was more investigated by Zhang in [34J . Since in this procedure the factorization 
is performed in backward direction, we call it BFAPINV (Backward FAPINV) algorithm. In 
[33], Zhang proposed an alternative procedure to compute the factorization in the forward 
direction, which we call it FFAPINV (Forward FAPINV) algorithm. In |14j . Lee and Zhang 
have shown that the BFAPINV algorithm is well-defined for M-matrices. It can be easily 
seen that this is correct for the FFAPINV algorithm as well. In [31], Salkuyeh showed that 
the FFAPINV algorithm with a simple revision may be used for the nonsymmetric positive 
definite matrices, free from breakdown. 

In this paper, we show that the FFAPINV algorithm is free from breakdown for H- 
matrices and we propose a technique for computing an ILU preconditioner based on the 
FAPINV algorithm at no extra cost. 

This paper is organized as follows. In section 2, we review the FFAPINV algorithm. 
Section 3 is devoted to the main results. Numerical experiments are given in section 4. 
Finally, we give some concluding remarks in section 5. 



2. A review of the FFAPINV algorithm 

Let W and Z T be lower unitriangular matrices and D be a diagonal matrix. Also, 
suppose that 

WAZ = D~ 1 . (2) 

In this case, we term W, Z and D, the inverse factors of A = (ctij). Consider W = 
(wj , wj, ■ ■ ■ , Wn) T , Z = (z\, Z2, ■ ■ ■ , z n ) and D = diag(di, &%, . . . , d n ), in which w^s and 
ZiS are the rows and columns of W and Z, respectively. Using Eq. ([2]) we obtain 

mA Zj = \ f i= J: (3) 

From the structure of the matrices W and Z, we have 

3-1 



z 1 = e 1 , Zj = ej -^af'zi, j = 2,...,n, (4) 
i=l 



3-1 

' ' '' ^ ; ^"'<- i 2 "• ( 5 ) 



i=l 



for some a\ ' s and p\ s, where ej is the jth column of the identity matrix. 
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First of all, from ([3]) we see that 

* 1 



z[Azi ejAei a u 

Now let 2 < j < n be fixed. Then from Eqs. ([3]) and @ and for k = 1, . . . , j — 1, we have 

= w k Azj = w k Aej - y^^aV'wkAzi = w k A*j - a> k w k Az k = w k A^ - a k ° ^ — , 

i=l fc 

where A*j is the jth column of A Therefore 



In the same manner 



otf' = diWiA*j, i = l,...,j-l. 



0% — diAj^Zi, i — 1, 



where A,-* is the jth row of A. Pre-multiplying both sides of (|3J) by WjA yields 

1 



(1: 



w j A*j 

Putting these results together gives the following algorithm (FFINV for Forward Factored 
INVerse) for computing the inverse factors of A. 

Algorithm 1: FFINV algorithm (vector form) 

1. z\ := e\, wi := ej and di := A- 

2. For j = 2, . . . , n, Do 

o _ T 

6. Zj :— ej\ Wj :— 

4. For i = 1, . . . , j — 1, Do 

5. a;p := diWiA^f, f3f' := d+Aj+Zi 

(i) (i) 

6. Zj := Zj - af'zi; Wj := Wj - /3f'Wi 

7. EndDo 
9. EndDo 

10. Return W = [wj , • • • , w^] T , Z = [zi, • • • , z n ] and D = diag(di, ■ ■ ■ ,d n ). 

Algorithm 1, is the vector form of the FFINV algorithm. It can be easily verified that 
this algorithm is equivalent to the following algorithm (see [34] ) . Moreover, the values of 

(i) (i) 

a/ 's and /3l JJ, s are the same in both Algorithms 1 and 2. In this algorithm, we assume 
Wj = (vjji,Wj 2 , . . • , w jn ) and Zj = (zij,z 2 j, z nj ) T . 

Algorithm 2: FFINV algorithm 

1. W := I-nxm Z '.= I n xn- 

2. For j = 1, . . . , n, Do 



3. 


For i = j - 


1,...,1, Do 


4. 


#> = ■ 


, 1 
^ji ' 2—/k=l ^jk^ki 


5. 


& = 




6. 


EndDo 




7. 


For i = 1, . 


Do 


8. 






9. 


EndDo, 




10. 


dj = 1/ (ajj 


+ Efc-1 w jk a kj) 


11. 


For i = j - 


1,...,1, Do 


12. 


»?» = 




13. 


<*« = 




14. 


EndDo 




15. 


For i = 1, . 


Do 


16. 


Z ij = ' 




17. 


EndDo 




18. 


EndDo 





19. Return = (ifjj), ^ = (%) and Z) = diag(d\, • • • , <i n ). 

This algorithm computes the inverse factors W, Z and D such that Eq. ([2]) holds. 
Therefore, we have A^ 1 = ZDW. This shows that the inverse of A in the factored form 
can be computed by this algorithm. A sparse approximate inverse of A is computed by 
inserting some dropping strategies in Algorithm 2. A dropping strategy can be used as 
follows. The dropping strategy is applied in four places in the algorithm. In step 5, if 

| Pi \ < t, then 1$?' := 0, and in step 8, if \wji\ < r, then Wji := 0. In the same way 

(i) 

a>i in step 13 and in step 16 are dropped when their absolute values are less than 
tolerance r. Hereafter, FFAPINV algorithm refers to FFINV algorithm with this type of 
dropping. Here we mention that the dropping strategy proposed in [33] is slightly different 
from our dropping strategy. In the next section, we show that the FFAPINV algorithm is 
well-defined for H-matrices. 



3. Existence of FFAPINV algorithm for H-matrices 

First we state the following theorem. 

Theorem 1. Assume that A is an M-matrix. Let W and Z be the inverse factors of A 
computed by FFINV algorithm, i.e., WAZ = D^ 1 . Also suppose that and W and Z be the 
inverse factors of A computed by the FFAPINV algorithm, i.e., WAZ ~ D . Then 

w >W > 0, Z>Z>0, 
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1 1 

-r>^>0, j = l,2,...,n, 



where D = diag(di, . . . , d n ). 



Proof. The proof of this theorem is quite similar to Proposition 2.2 in [14] and is omitted 
here. □ 

This theorem shows that the FFAPINV algorithm is well-defined for M-Matrices. We 
mention that this is correct for the BFAPINV algorithm, as well (see |14j). 

Theorem 2. Let A be an H-matrix and A be its comparison matrix. Let also A^ 1 = ZDW 
and Ar x = ZDW be the computed inverse in the factored form by FFINV algorithm for A 
and A, respectively. Then 

,1, 1 

di di 

Proof. The elements Wji and z^ may be assumed as rational functions 

W jl = Fjl( a U: ■ ■ ■ > a j-l,n, di, . . . , dj-i), 

z ij = Gij(aii,...,a n j-i,di,...,dj-i). 

In fact, Wji (zij) is a rational function of first j — 1 columns (j — 1 rows) of A and first j — 1 
diagonal entries of D. In the same way, Wji {zu) is a rational function Fji (Gy) of first j — 1 
columns (j — 1 rows) of A and first j — 1 diagonal entries of D. Let us also assume that 

Wji = Fji(an, . . . , %-i )T i) \di\, ■ ■ ■ , 
% = Gij{d\i, . . . ,a„j_i, . . . , 

This means that Wji and zij are computed similar to Wji and %, with . . . , \dj-i \ instead 
of d±, . . . , dj—i. By induction, we prove that 

b) w kt > w k t > 0, 

c) f t k) <l[ k) <0, (6) 

d) hk > ztk > 0, 

e) 4 k) <u { t k) <0, 

for k = 1, ... ,n and t < k — 1. Note that l^ k \ & , ui k ' and u[ are defined similar to 
Wji and Wji. For fc = 1, there is nothing to prove. Now, let all of these relations hold for 
every k < j — 1. We show that all of them are correct for k = j, as well. From step 4 of 
Algorithm 2, for every t < j — 1, we have 

t-l 

Z/^ = cijt + ajiZit- 
i=l 
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From the hypothesis for every t < j — 1, we have za > % > 0. Therefore 

f ) = %t + E*=i %^it < %t + E*=i = ? J) < o- ( 7 ) 

Also, from steps 7-9 of Algorithm 2, we have 

^i* = -it ] dt - ^2 Vf'diW iu t = 1, . . . , j - 1. 

i=t+l 

From (J7|) and the hypothesis, for every k < j — 1, we have 

4 > 141 > o -z£°4 > 

w)fc t > *h > -l^d k w k t> -l^\dk\w k t, t < k. 
Hence, we conclude that 

i-l j-i 

i=t+l i=t+l 

In the same way, it can be verified that 

i-l i-l 

■u^ = d t j + ^2 dijWu < d t j + ^2 a ij™tj = < 0, 
i=l i=i 

i-i i-i 

% = -u[ j) d t - ^2 ^PdiZu > -u[ j) \d t \ - ^2 \di\z u = %jt > 0. 

i=t+l i=t+l 

Now, we consider two cases ajj > and ajj < 0. If ajj > 0, then from Wjt > Wjt > and 
dtj < 0, we conclude that 

x i-i i-i 

J = %i + X] - a w + atjWjt- 

d j t=i t=i 

Consider the polynomials of Wjt and Wjt- It is easy to see that the corresponding terms of 
these polynomials have the same absolute values. In other words, they may differ only by 
the sign. On the other hand, every term of Yli=l dtjWjt, when considered as a polynomial 
in elements of A, is nonpositive, since all terms of w are nonnegative. Therefore, it is less 
than or equal to J2t=i a tjWjt- Since, it is enough that one of its terms to be nonnegative. 
Putting these results together indicates that 

1 t$ 1 

y < ajj + ^2 at3™it < ajj + ^2 a ti w 3t = ~T ■ (8) 
d i t=i t=l d i 

If ajj < 0, then from Wjt > Wjt > and dtj < 0, we have 

1 i-i i-i 

d i t=i t=i 
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Similar to the case ajj > 0, one can conclude that 

j-i i-i i- 1 



J J J 

y > -ajj - atjWjt > -a# + ^ ayiity = a# + ^ a ti m, t = — . (9) 
d j t=i t=i t=i a J' 



Now from Eqs. ([8]) and (jSJ), we have 

11 

dj ~ d/ 

This together with Theorem 1 give the desired result. □ 

Theorem 3. Let A be an H-matrix and A be its comparison matrix. Let also A~ x ~ ZDW 
and Ar x ~ ZDW be the factored approximate inverse computed by FFAPINV algorithm 
for A and A, respectively. Then 

|i|>i>o. 

di d t 

Proof. The proof is quite similar to the proof of the previous theorem. Indeed, the proof 
can be done by induction and noting that all of the inequalities have been obtained by 
comparing both sides of the inequalities term-by-term. □ 

From this theorem, the following corollary can be easily concluded. 

Corollary 1. Let A be an H-matrix and D = diag(di, . . . , d n ) be the diagonal matrix 
computed by Algorithm 1. Then, the sign of ajj and dj are the same. 



3. An ILU preconditioner based on the FFAPINV algorithm 

We first state and prove the following theorem. 

Theorem 1. Let W and Z be the inverse factors of A in Eq. H)). Also suppose that 
L := W~ x and U := Z^ 1 . Therefore A = LD~ l U is the LDU factorization of A and for 
i < 3 

(i) (i) 

in which a\ and /?• are computed in steps 5 and 13 of Algorithm 2. 

Proof. Let WAZ = D~ l . Then A = W~ X D~ X Z~ X is the LDU factorization of A. From 
= U we conclude ZU = L, where / is the identity matrix. By equating the jth column 
of both sides of ZU = I, we deduce that 

i-i 

Zj = ej - ^UijZi, (10) 
i=i 
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where ej is the jth column of /. Since Zi, i = 1, 2, . . . ,j — 1, and ej are linearly independent 
vectors, then relation (jlOp together with results in = U%j. In the similar way we 
can prove that Lji = j3\ f . □ 

By the above discussion we propose the next algorithm that computes an ILU factor- 
ization of A as by-product of FFAPINV process. We term this algorithm ILUFF (refer to 
an ILU preconditioning based on FFAPINV algorithm). The algorithm is as follows: 

Algorithm 3: ILUFF algorithm 

1. Set L = U = I nxn , Z\ := e\, wi := ej and d± := ^_ 

2. For j = 2, . . . , n, Do 



o _ _ T 

6. zj .— ej, Wj .— ej 

4. For i = 1, . . . , j — 1, Do 

5 . \J I j . — di TVi ./I % j 

6. If \Uij\ > t, then Zj := Zj — UijZ; L 

7. Drop entries of zj whose absolute values are smaller than r 

8. EndDo 

9. For i = 1, Do 

1 . Ij j >i • — d-i -A-j ^ %i 

11. If \Lji\ > t, then Wj := Wj — LjiWi 

12. Drop entries of Wj whose absolute values are smaller than r 

13. EndDo 

15. EndDo 



16. Return L = (Lji), U = {Uij) and D = diag(di, d 2 , ■ • • , d n ) (A « LD^U) 

Some consideration can be given here. Obviously, when the matrix A is symmetric then 
W = Z T and the computations will be halved. In the case that the matrix A is symmetric 
positive definite then we have 

^ = ^- >0 ' 

and the algorithm is free from breakdown. This is true for nonsymmetric positive defi- 
nite matrices, as well [31]. Therefore, if the matrix A is positive definite (symmetric or 
nonsymmetric) then we can use 

' / ' : ?k' 

in step 14 of Algorithm 3. 
4. Numerical experiments 
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In this section, we have used the ILUFF as the right preconditioner to solve the linear 
system of equations (H|) with GMRES(50) method. The ILUFF code is written in Fortran 
77. But the GMRES(50) code in Sparskit package [30] has been used. In the first part of 
the numerical experiments we used 38 nonsymmetric test matrices. All of these matrices 
have been taken from the University of Florida Sparse Matrix Collection [U] and none of 
them are positive definite. In all the experiments whenever a zero pivot has occurred, then 
we have replaced the zero by the square root of the machine precision. The machine, we 
used for the experiments, has one quad Intel(R) CPU and 8GB of RAM memory. The 
initial guess for the iterative solver was always a zero vector. The stopping criterion used 
was 

Kl| 2 < 10 -i 0; (n) 



F0II2 

where is the residual of the unpreconditioned system in the kth iterate. We have consid- 
ered the exact solution as the vector e = (1, • • • , 1) T and the vector b = Ae. We have used 
the Multilevel Nested Dissection reordering as a preprocessing |18j . To implement ILUFF 
algorithm, it is clear that matrix A should be accessed row-wise and column-wise. In our 
implementations, we have used just the Compressed Sparse Row (CSR) format of storage 
[29j to traverse A row-wise. For the column-wise traverse of A, we have used the linked 
lists pi]. 

In Table [H properties of test matrices and the convergence results of the GMRES(50) 
method without preconditioning have been reported. In this table, n and nnz are the 
dimension and the number of nonzero entries of the matrix, respectively. Its stands for the 
number of iterations and Time is the iteration time which has been computed by dtime 
command. The times are in second. In this table, a + symbol means that the stopping 
criterion has not been satisfied after 10,000 number of iterations. 

Table El includes properties of ILUFF preconditioner and results of right preconditioned 
GMRES(50). Considering Algorithm 3, r is the tolerance parameter to drop entries of L, U 
and W, Z factors. In this table, the parameter r has been set to 0.1 for all the test matrices. 
Suppose that nnz(A), nnz(L) and nnz(U) be the number of nonzero entries of matrix A and 
L and U factors of ILUFF preconditioner. In the implementation of ILUFF preconditioner, 
we have merged the D factor into the U factor. Therefore, parameter density in Table [21 
is defined as: 

nnz(L) + nnz(U) 

density = -—- . 

nnz(A) 

In this table, Ptime stands for the preconditioning time and It — Time is the iteration 
time to solve the preconditioned system. Ttime is defined as the sum of Ptime and It — 
Time. In this table, Its is the number of iterations of the GMRES(50) to solve the right 
preconditioned system. Numerical results presented in Table 2, show that the proposed 
preconditioner greatly reduces the time and iterations for convergence. 

For the second part of the numerical experiments we consider the matrix atmosmodj . 
This matrix belongs to the Group Bourchtein of matrix collection [9]. Dimension of this 
matrix is 1,270,432 and it has 8,814,880 number of nonzero entries. Consider system 
([1]) when the coefficient matrix is atmosmodj . Also suppose that b = Ae in which e = 
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Table 1: Test matrices properties together with results of GMRES(50) without precondi- 
tioning. 



1 rAiir* / \ .1 ott*iv 

Vjioup / IVldll IX 


71 


TITIZ 


Time 
_L 1111c 


lib 


HB/is_ 183_ 1 


183 


998 


0.00 


38 


HB/ts_ 183_ 6 


183 


1000 


0.01 


36 


Simon/racfskyl 


3242 


293409 


2.76 


3588 


Simon/racfsky2 


3242 


293551 


3.65 


4790 


Simon/racfsky5 


6316 


167178 


0.24 


306 


Simon/ ractskyb 


3402 


130371 


0.68 


1385 


TV T * i / /"I 11 1 O 

Muitc/Chebyshev3 


4101 


36879 


+ 


+ 


(J berwoliacn/nownietero 


9669 


67391 


+ 


+ 


Rajat/rajat03 


7602 


32653 


+ 


+ 


HB / shcrman3 


5005 


20033 


+ 


+ 


TT / 1 

Hamm / mcmplus 


17758 


99147 


7.22 


3878 


FEMLAB / poisson3Da 


13514 


352762 


0.86 


342 


Botonakis/FEM_ 3D_ thcrmall 


17880 


430740 


0.88 


283 


FEMLAB / poisson3Db 


13514 


352762 


12.98 


620 


/~\1 11" 1/1' 1 r\ 

(J bcrwolf ach/chipcoolO 


20082 


281150 


+ 


+ 


/~\1 li" 1/1" 11 

(J bcrwoliacn/chipcooll 


20082 


281150 


+ 


+ 


A / 1 1 

Averous/epbl 


14734 


95053 


2.08 


1432 


A / 1 O 

Averous/epbz 


25228 


175028 


2.81 


908 


Wang/wang3 


26064 


177168 


2.34 


803 


Wang/wang4 


26068 


177196 


+ 


+ 


IBM_ Austin/coupled 


11341 


97193 


+ 


+ 


Simon/venkatOl 


62424 


1717792 


+ 


+ 


Sandia/ASIC_ lOOks 


99190 


578890 


23.85 


1887 


Hamm/hcircuit 


105676 


513072 


+ 


+ 


Norris/lung2 


109460 


492564 


+ 


+ 


IBM_ EDA/dcl 


116835 


766396 


+ 


+ 


IBM_ EDA/dc2 


116835 


766396 


+ 


+ 


IBM_ EDA/dc3 


116835 


766396 


+ 


+ 


IBM_ EDA/trans4 


116835 


749800 


+ 


+ 


IBM_ EDA/trans5 


116835 


749800 


+ 


+ 


Botonakis/FEM_ 3D_ thcrmal2 


147900 


3489300 


18.32 


652 


QLi/crashbasis 


160000 


1750416 


10.09 


437 


FEMLAB /stomach 


213360 


3021648 


11.64 


344 


Sandia/ASIC_ 320ks 


321671 


1316085 


10.39 


201 


Sandia/ASIC_ 680ks 


682712 


1693767 


13.72 


84 


Bourchtein/atmosmodd 


1270432 


8814880 


241.92 


707 


Bourchtcin/atmosmodl 


1489752 


10319760 


166.14 


415 
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Table 2: Properties of ILUFF preconditioner with r = 0.1 and right preconditioned GMRES(50) 
method 



Matrix 


density 


Ptime 


It-Time 


Ttime 


Its 


fs_ 183_ 1 


0.55 


0.81 


0.00 


0.81 


10 


fs_ 183_ 6 


0.54 


0.73 


0.00 


0.73 


10 


raefskyl 


0.05 


0.75 


1.03 


1.78 


1092 


raefsky2 


0.08 


0.75 


1.45 


2.20 


1419 


raefsky5 


0.20 


0.75 


0.01 


0.76 


11 


raefsky6 


0.15 


0.75 


0.01 


0.76 


12 


Chebyshev3 


0.33 


0.75 


0.09 


0.84 


245 


flowmctcr5 


0.74 


0.76 


2.37 


3.13 


2106 


rajat03 


0.78 


0.82 


0.25 


1.07 


390 


shcrman3 


0.83 


0.73 


0.82 


1.55 


1747 


mcmplus 


0.39 


0.77 


0.89 


1.66 


376 


poisson3Da 


0.18 


0.75 


0.55 


1.30 


180 


FEM_ 3D_ thermall 


0.22 


0.76 


0.23 


0.99 


53 


poisson3Db 


0.18 


1.06 


9.40 


10.46 


395 


chipcoolO 


0.35 


0.78 


1.22 


2.00 


321 


chipcooll 


0.35 


0.75 


1.70 


2.45 


446 


epbl 


0.78 


0.73 


1.11 


1.84 


547 


epb2 


0.57 


0.75 


0.83 


1.58 


209 


wang3 


0.85 


0.76 


0.84 


1.60 


201 


wang4 


0.55 


0.80 


1.13 


1.93 


279 


coupled 


0.48 


0.81 


0.21 


1.02 


149 


vcnkatOl 


0.34 


1.07 


1.67 


2.74 


90 


ASIC_ lOOks 


0.79 


1.01 


0.35 


1.36 


23 


hcircuit 


0.76 


0.99 


8.95 


9.94 


513 


lung2 


1.03 


1.00 


5.67 


6.67 


304 


del 


0.65 


35.26 


4.75 


40.01 


234 


dc2 


0.64 


34.06 


2.97 


37.03 


143 


dc3 


0.64 


33.75 


8.64 


42.39 


451 


trans4 


0.62 


21.95 


2.36 


24.31 


128 


trans5 


0.63 


21.80 


7.52 


29.32 


397 


FEM_ 3D_ thcrmal2 


0.22 


1.19 


1.41 


2.60 


41 


crashbasis 


0.58 


1.19 


1.82 


3.01 


59 


stomach 


0.22 


1.26 


2.79 


4.05 


72 


ASIC_ 320ks 


0.66 


1.45 


4.02 


5.47 


71 


ASIC_ 680ks 


0.61 


1.97 


0.68 


2.65 


6 


Bourchtein/atmosmodd 


0.64 


4.02 


200.66 


204.68 


503 


Bourchtcin/atmosmodl 


0.84 


5.14 


100.32 


105.46 


209 
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Table 3: Properties of ILUFF preconditioner and right preconditioned GMRES(50) for the matrix 
atmosmodj 



T 


density 


Ptime 


It-Time 


Ttime 


Its 


0.1 


0.64 


4.04 


240.81 


244.85 


603 


0.05 


0.76 


4.42 


232.66 


237.08 


561 


0.01 


1.04 


5.17 


198.06 


203.23 


464 



(1, • • • , 1) is the exact solution. We term this artificial system as the atmosmodj system. 
Without preconditioning, GMRES(50) method for this system converges in 1312 number 
of iterations in about 447.66 seconds. In Table El we have reported the results of ILUFF 
preconditioner and the right preconditioned GMRES(50) method for this system. In this 
table, density, Ptime, It — Time, Ttime and Its have the same meaning as in Table [2j All 
the Ttime and Its in this table are less than numbers 447.66 and 1312, respectively. 

In figure [H we have drawn four graphs related to atmosmodj system. In this figure, we 
take an in-depth look at the results of Table O The graph with solid line is devoted to the 
case that GMRES(50) method without preconditioning is applied to solve the atmosmodj 
system. This graph gives the logarithm of the fraction |l* E f° r each iterate of the GM- 
RES(50). The three other graphs, illustrated by the dashed, dotted and dashed-dotted 
lines, give the above mentioned logarithm for each iterate of the right preconditioned 
GMRES(50) method. For these three graphs, the right preconditioner is the ILUFF and has 
been computed with r equal to 0.1, 0.05 and 0.01. The shape of the graphs first indicates 
that the ILUFF preconditioner is a good tool to decrease the number of iterations of the 
GMRES(50) method, applied for atmosmodj system. It also shows the fact that the less 
the parameter r gives the better the quality of the ILUFF preconditioner for atmosmodj 
system. 



5. Conclusion 

In this paper, we have proposed an ILU preconditioner based on the Forward FAPINV 
algorithm say ILUFF which is free from breakdown for nonsymmetric positive definite and 
H-matrices. Numerical results presented in this paper show that the new preconditioner is 
very robust and effective. 
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